Structures in magnetohydrodynamic turbulence: detection and scaling 
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We present a systematic analysis of statistical properties of turbulent current and vorticity struc- 
tures at a given time using cluster analysis. The data steins from numerical simulations of decaying 
three-dimensional (3D) magnetohydrodynamic turbulence in the absence of an imposed uniform 
magnetic field; the magnetic Prandtl number is taken equal to unity, and we use a periodic box 
with grids of up to 1536'' points, and with Taylor Reynolds numbers up to 1100. The initial condi- 
tions are either an X-point configuration embedded in 3D, the so-called Orszag-Tang vortex, or an 
Arn'old-Beltrami-Childress configuration with a fully helical velocity and magnetic field. In each 
case two snapshots are analyzed, separated by one turn-over time, starting just after the peak of 
dissipation. We show that the algorithm is able to select a large number of structures (in excess of 
8, 000) for each snapshot and that the statistical properties of these clusters are remarkably similar 
for the two snapshots as well as for the two flows under study in terms of scaling laws for the 
cluster characteristics, with the structures in the vorticity and in the current behaving in the same 
way. We also study the effect of Reynolds number on cluster statistics, and we finally analyze the 
properties of these clusters in terms of their velocity-magnetic field correlation. Self-organized crit- 
icality features have been identified in the dissipative range of scales. A different scaling arises in 
the inertial range, which cannot be identified for the moment with a known self-organized criticality 
class consistent with MHD. We suggest that this range can be governed by turbulence dynamics 
as opposed to criticality, and propose an interpretation of intermittency in terms of propagation of 
local instabilities. 

PACS numbers: 47.65.-d,47.27.Jv,96.50.Tf,95.30.Qd 



I. INTRODUCTION 



As the resolving power of experimental instrumenta- 
tion increases, turbulent flows as they occur in geophysics 
and astrophysics are being examined with more accu- 
racy, and the multiplicity of scales in interactions be- 
comes more apparent. To take an example, the modal 
distribution of energy in the Solar Wind has been known 
for a long time to follow a power law close to the Kol- 
mogorov prediction for incompressible fluid turbulence 
[ij (see [31 for reviews), although its physical environment 
is infinitely more complex than what was first envisaged 
by Kolmogorov, involving magnetic fields and coupling 
to acoustic and whistler modes, to name but a few phe- 
nomena at play. There are also numerous observations 
of spatially correlated turbulent structures and flows in 
Earth's magnetosphere. Recent observations from Clus- 
ter and THEMIS multi-spacecraft missions 0, |j| provide 
a sophisticated physical picture of a variety of significant 
effects, for example, intermittent (spatially sparse) struc- 
tures and transient plasma transport associated with re- 
connection in the tail plasma sheet and at the dayside 
magnetopause, formation of shocks and small-scale mag- 
netic filaments, Kelvin-Helmoltz vortices and coherent 
structures viewed as Alfvenic turbulence Q, as well as 
other effects. Signatures of MHD turbulence are also 
found, e.g., in the magnetosphere of Jupiter [6|, and in 
the interstellar medium [3]. 



Of the many features now being resolved in the ob- 
servations, intermittency in turbulence is of critical im- 
portance as it is linked with magnetic energy conver- 
sion and dissipation in solar-terrestrial plasmas. Two 
well-known examples of this link are fiaring activity in 
the solar corona, and magnetospheric substorms in the 
tail plasma sheet of Earth's magnetosphere. In both 
cases, free magnetic energy is released through spatially 
localized reconfiguration of the plasma geometry which 
is significantly affected by MHD turbulence. Intermit- 
tent magnetic structures in the solar corona can gener- 
ate multiple tangential discontinuities leading to energy 
avalanches and strongly inhomogeneous dissipation [8|- 
[lOj . An enhanced intermittency often reflects the forma- 
tion of an unstable magnetic topology. The latter has 
been explored in detail when examining data from solar 
active regions which reveal precursory dynamics of inter- 
mittent measures prior to large solar flares [ll|. MHD 
intermittency is also likely to be a major factor defin- 
ing initial locations of magnetic reconnection events in 
the nearly coUisionless plasma of the Earth's magneto- 
tail |12l - [l4| . It can be a triggering mechanism for a va- 
riety of instabilities of p lasma behavior at both kinetic 
and MHD scales [IJ, UM responsible for multiscale parti- 
clc precipitation in the night-time auroral region [l^, |l3l . 
The timing; position, and energy output of such events 
- as well as the structure of the Solar Wind mediating 
their interaction - are largely unpredictable, refiecting 



the stochastic nature of the underlying fluid dynamics. 

The intermittent structures associated with dissipation 
in turbulent flows are difficult to detect because they re- 
side mostly at small scale, in thin current and vorticity 
sheets, and are entangled with ambient plasma flows in a 
topologically complex way. These structures are dynam- 
ically important as they participate in the local heating 
of the medium, for example through localized magnetic 
reconnection events, and are important not only in as- 
trophysics but also in laboratory plasmas (see, e.g., [18|). 
Intense and localized dissipative structures in MHD flows 
have also been obtained numerically, in both two (2D) 
[19| and three space dimensions (3D) [20|. Intermittent 
dissipation in MHD simulations is found to be typically 
stronger than that in neutral fluids (i.e., with a stronger 
departure from self-similarity), and it can vary in time, 
as observed in solar active regions [ll|, [23 . A striking fea- 
ture of these structures, as found in many simulations, is 
a high degree of correlation between the velocity and the 
magnetic field both globally [21| and locally ^^. How- 
ever, the precise relationship between these structures 
and the global statistical properties of the flow is not 
well understood, although it is known that in 2D and 
in 3D the structures can interact with the underlying 
turbulent flow to affect properties such as the global dis- 
sipation through local processes like reconnection. Fur- 
thermore, at high Reynolds numbers these structures can 
have complex geometries, e.g., roll- up and fold as ob- 
served in the Solar Wind [29[ and in direct numerical sim- 
ulations (DNS) [30], complicating significantly attempts 
to make a connection between structures and statistical 
properties. 

Ensemble-based description of the geometry of inter- 
mittent dissipation is an important issue as a turbulent 
flow is characterized not only by the structures that de- 
velop within, but also by the statistical properties of the 
flow as a whole. For fluids in 2D, the statistics of vor- 
tices have been studied in detail 12J| and a relationship 
between vorticity and stream function has been found 
which can be ascribed to a distribution of signed vor- 
tices l25i using a maximum entropy principle [26| (see 
also [231 for a recent analysis). In three dimensions, the 
situation is much more complex but in one specific case, 
a Kolmogorov spectrum for the energy has been obtained 
analyticall y fr om the dynamics of the stretching of a spi- 
ral vortex |28j . However, while high- intensity dissipative 
structures in 3D MHD have been successfully studied 
for a number of years [3l|, mostly through threshold- 
ing and visualization of current and vorticity (see, e.g., 
[33-[35| in 3D, and [3a] for a thorough study of recon- 
nection events in 2D), we are not aware of any ensemble- 
based studies of turbulent structures observed at inter- 
mediate to small intensity levels in high Reynolds num- 
ber 3D MHD. A quantitative analysis, which requires 
the development of new software tools, is particularly 
important in the wake of two overarching developments: 
the emergence of petascale computers that will produce 
vast amount of data and detailed point- wise information 



about the relevant dynamical variables and their deriva- 
tives, as well as planned in situ observations, in particular 
those in association with the upcoming NASA's Magne- 
tospheric Multiscale (MMS) mission, which will investi- 
gate the role of turbulence and other cross-scale phenom- 
ena in fast magnetic reconnection. 

Therefore, in this paper we propose a new methodology 
for analyzing cross-scale behavior of three-dimensional 
MHD turbulence enabling the detection of multiple dis- 
sipative structures at arbitrary intensity levels. We use 
our tools to extract current and vorticity structures in 
numerical simulation outputs with two distinct types of 
initial conditions. The results obtained show the exis- 
tence of robust scaling behavior in both the inertial and 
dissipative regimes of scales in the turbulent fluids we 
study. The reported scaling exponents shed new light 
on the role of the initial conditions and of the Reynolds 
number in the formation of intermittent dissipative struc- 
tures in current and vorticity flelds. Finally, our analysis 
supports the possibility of self-organized critical behavior 
for some of the small-scale structures we detect with our 
algorithm. 



II. METHODOLOGY 

A. Equations and flows 

The incompressible MHD equations in dimensionless 
Alfvenic units and in the absence of forcing read: 



c'tv + V • Vv = -/9(7^ VT' + j X b + i/V^v, 


(1) 


a^b = V X (v X b) + jyV^b , 


(2) 



with V, b being the velocity and the magnetic flelds, 
j = V X b the current density, V the pressure, po=l 
the constant density, and V ■ v = V • b ~ 0. When 
the viscosity v and the magnetic resistivity 77 are both 
equal to zero, the energy Et = (v^ + b'^) /2, cross hclic- 
ity He = (v • b) /2, and magnetic helicity Hm = (A • b) 
(where A is the vector potential, b = V x A) are con- 
served. Equations dl])-© have been solved in a three- 
dimensional box with periodic boundary conditions and 
a pseudospectral code dealiased by the 2/3 rule; /cmin = 1 
for a box of length Lq = 27r, and with N regularly 
spaced grid points, this leads to a maximum wavenum- 
ber fcmax = iV/3. At all times, fcu/femax < 1, with ko 
denoting the dissipation wavenumber, in order to ensure 
an accurate numerical computation down to the smallest 
resolved scale. 

Two types of flows are studied in this paper, which 
were computed on regular grids ranging from 512^ to 
1536'^ points (see Table U for a brief presentation of the 
runs; see also [33, [S^l where these runs are described in 
the context of a study of the general properties of MHD 
turbulence) . In the first flow (runs I and II) , the flelds are 
constructed from a superposition of Arn'old-Beltrami- 
Childress (ABC) flows (see, e.g., [331), ^^ wavenumbers 



TABLE I: Nomenclature of runs with A^ the Unear grid 
resolution, and v and r\ the viscosity and magnetic diffu- 
sivity, respectively. The Reynolds number Kv = UoL{)/u 
and Taylor Reynolds number R\ — UqX/u (with Uo the 
r.m.s. velocity, Lo a characteristic large-scale and A — 
2it[Et / j k^ ET{k)dkY'^ the Taylor scale based on the total 
energy spectrum) are both evaluated at peak of dissipation. 



with ath G [Ij 3]: 



Run 


Type 


N 


V = rj 


Rv 


Rx 


I 

II 
III 


ABC 

ABC 

OT 


512 
1536 
512 


6. X 10-* 
2. X 10"* 
7.5 X 10"'' 


3100 
9200 
3300 


630 
1100 
880 



1 to fc 



tuations are added with a spectrum k ^ exp[— 2(fc/fco)]^ 



to which smaller-scale random fluc- 
'exp[ 
for fc > 3 (see |37|). The phases of the modes with 
A; > 3 arc chosen from a Gaussian random number gen- 
erator in such a way that the initial cross-correlation of 
the two fields is negligible: initially, Ey = Em = 0.5, 
He « 10-^ and Hm = 0.45. 

Another flow we compute (run III) is that of the so- 
called Orszag-Tang vortex (OT hereafter) generalized to 
MHD [32[ (sec also [sj] and references therein) . This flow 
has been studied at length in two space dimensions for its 
reconnection properties (it has a magnetic X-point cen- 
tered at a stagnation point of the velocity) ; its general- 
ization to 3D is straightforward, with a simple sinusoidal 
variation in the z direction. Initially, Ey = Em = 2, 
He = 0.41 and Hm = 0. Note that the two types of 
initial conditions differ in their invariants: the OT flow 
has zero magnetic hclicity and a sizable cross correlation, 
whereas for the ABC flow, it is the opposite. 



B. Cluster detection 

Turbulent flows exhibit small-scale structures with 
strong gradients in the vicinity of which dissipation takes 
place. In principle, detection of structures can be done on 
any physical variables but more essential to a turbulent 
flow with its complex small-scale behavior are vorticity 
and current. Indeed, the primary channels of spatially 
localized energy dissipation in a resistive MHD fluid are 
the Joule heating, proportional to the squared current 
density p = |V x b|^, and the kinetic dissipation that 
can be characterized indirectly by the squared vorticity 
uP' = |V X v|^ (note that the local dissipation of kinetic 
energy is proportional to the symmetric part of the ve- 
locity gradient matrix, the difference stemming from the 
fact that V is a vector whereas b is an axial vector). 

Our analysis is thus focused on 3D arrays containing 
the values of j"^ and w^ for each of the A^'^ grid nodes of 
the simulations listed in TablelH A grid node is treated as 
belonging to a dissipativc structure if the dissipation in 
this node, expressed in terms of either field, exceeds the 
level of ath standard deviations above the mean value, 



,■2 

Jth 



'^th 



{f)+atH^{j^)-{{P))\ (3) 

(..2)+a„V(^')"((^'))'. (4) 

and where (•) denotes averaging over the entire simu- 
lation volume as before. Intermittent dissipation struc- 
tures in the p field (current sheets) are defined as spa- 
tially connected sets of grid nodes satisfying the condi- 
tion p > j^^\ the structures in the w^ field (vorticity 
sheets) are defined similarly, based on w^ > w^^. Our 
goal is to separate these structures from the background 
and to study their individual properties such as linear 
size, volume, or internal dissipation rate, for subsequent 
ensemble-based surveys. 

In order to overcome inherent memory limitations of 
standard cluster analysis algorithms such as, e.g., the 
Label Region function of the Interactive Data Language 
(IDL), we have developed a new technique enabling fast 
decomposition of multidimensional data arrays into sets 
of dissipative structures while dramatically reducing the 
amount of stored information. 

The first step of our technique consists of building a ta- 
ble of contiguous intervals (activation sites) along a cho- 
sen direction ("scanning direction") where the studied 
data field exceeds the detection threshold. We tabulate 
boundaries of such intervals for all relevant positions in 
the d — 1 coordinate space orthogonal to the scanning 
direction, d being the dimensionality of the original data 
set. By design, the choice of the scanning direction is 
arbitrary and does not affect the detection results. The 
second step is to find and label spatially connected clus- 
ters of activations using the "breadth-first search" prin- 
ciple to avoid backtracking of search trees representing 
individual clusters. We find that it is important to con- 
sider all of the 3^^ — 1 nearest neighbors of each grid node, 
including the diagonal neighbors, when identifying con- 
nected activations. Finally, the activation table is sorted 
a second time according to the cluster labels, in order to 
provide faster access to the detected structures. 

The output data array preserves complete information 
on the location and shape of all the contiguous regions in 
the simulation volume where the threshold condition is 
fulfilled. The size of this array is typically smaller by two 
order of magnitudes (depending on the threshold), com- 
pared to the original data. For the purpose of this paper, 
the described technique is used to identify structures in 
static snapshots of a turbulent fluid. However, it can 
be easily extended to higher-dimensional data sets repre- 
senting spatiotemporal dynamics of turbulent structures, 
with the time axis being a natural scanning direction. 



C. Data analysis tools 

The data sets analyzed in this paper are three- 
dimensional cubes {d ~ 3) described by x, y and z coor- 
dinates. The scanning direction is chosen parallel to the 




FIG. 1: [Color online) Top: Global view of dissipative clus- 
ters in the j^ field selected by our algorithm in the the ABC 
run I. The largest cluster has been removed to let see the 
intermediate size clusters; only one tenth of the remaining 
clusters are shown. Bottom: Zoomed-in view of two selected 
current clusters, showing strong curvature of the sheets; the 
vorticity behaves similarly. The lower right edge (direction 
marked with red in the color version) is parallel to the z axis 
chosen as the scanning direction. 




FIG. 2: [Color online) Two-dimensional visualization of the 
current at fifteen regularly spaced slices in the z direction with 
z £ [0,70] (in pixel units). A dissipative cluster made up of 
two twisted current sheets merging is shown in the middle of 
the box (highlighted in red in the color version). 



z axis. The activation tables represent lists of z inter- 
vals with [j(x, 2/, z)Y > Jth O'^ [^{^j Vi ^)Y > ^"iin ordered 
according to their projected position onto the x — y plane. 



The following primary parameters are computed and 
stored for each of the detected current and vorticity struc- 



turcs: 
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(10) 



Here, A^ is the set of all the grid node indices belonging 
to the ith dissipative structure, the subindices k and I 
label individual grid nodes, M{k) is the full set of the 26 
nearest neighbors of the fcth node, rfc is the dimension- 
less position vector of the fcth node, and finally 6 = 2tt/N 
is the grid spacing, uniform and isotropic. The primary 
parameters of the clusters are then as follows. Li is the 
linear size of the structure defined as the largest pair- 
wise distance (diameter) between the points tabulated in 
^iy L^^yz}i are the maximum dimensions of the struc- 
ture projected onto each spatial direction, Ri is the char- 
acteristic linear scale of the smallest Euclidean volume 
embedding the whole structure, Vi is the physical vol- 
ume occupied by the structure, and Ai is the area of 
its outer surface. Py^yi is the volume-integrated con- 
tribution of the structure to the kinetic and magnetic 
enstrophies, with the global corresponding quantities de- 
fined as ilj — J j^dV and 51^ — J uj^dV respectively, and 
where the integral is over the entire box. As mentioned 
before, except for a constant of proportionality {v = rf) 
these quantities are also the dissipation rates, and in the 
following we refer to P{j^uj}i as the "kinetic and magnetic 
dissipation" of the structure. 

In addition to these primary parameters, the following 
measures characterizing the geometry of the structures 
were used: 






v(rfc) ■b(rfc) 
|v(r,.)||b(rfc)|/,g^_ 
V^|[A,|2) 



(11) 

(12) 
(13) 



where Xi is the cosine of the local alignment angle be- 
tween the velocity and the magnetic field vectors aver- 
aged over all the positions within the ith structure. Hi 
and Ci are respectively the characteristic thickness of the 
structure and its topological complexity, both computed 
under the assumption of a sheet-like geometry which we 
validate later in the text. C is defined as the ratio be- 
tween the area of a circle with the diameter equal to the 



linear size Li of the structure, and the actual area of one 
of the sides of the structure. As follows from the def- 
inition, C increases if the structure has holes or other 
irregularities reducing A for a given i, and decreases if 
the structure has a curved shape ensuring a more efficient 
spatial filling for a fixed linear scale. The relative con- 
tribution of the second effect is expected to grow with 
L (larger structures tend to roll-up and fold more fre- 
quently than the smaller ones), while the first effect is 
nearly scale-invariant as we show below. 

Two major groups of statistics are invoked to quan- 
tify the scaling behavior of the detected structures. The 
first group includes a set of regression plots character- 
izing the geometric scaling of the structure parameters 
with respect to the linear size _L; the second group is 
represented by a set of probability distribution functions 
of structure parameters. Both types of statistics are ap- 
proximated by power-law dependencies: 



X{L) ex i^^. 



(14) 
(15) 



in which X denotes any of the parameters defined in Eqs. 
([5|)- ([T^ . p{X) is an estimated probability density func- 
tion of X, and Dx and tx are the geometric and the 
distribution scaling exponents, respectively; the latter 
are evaluated separately for the inertial and dissipative 
subranges of the fiows we study. The ranges of linear 
scales corresponding to these subranges are determined 
based on the behavior of the energy spectra. For runs 
I and III, the inertial range scaling is observed for wave 
numbers k E [5,30], an interval which corresponds to 
L £ [0.21, 1.3]. For run II, the inertial behavior is realized 
within the interval k G [5, 50] yielding L G [0.13, 1.3]. For 
the dissipative (sub-inertial) scaling regime, we choose 
L e [0.025, 0.18] in runs I and III and L e [0.012, 0.11] in 
run II. 

The inertial and dissipative scaling ranges as specified 
in terms of L are first applied to compute power-law fits 
describing the X{L) statistics. Next, the fits are used to 
evaluate the inertial and the dissipative scaling ranges of 
the remaining parameters, and to estimate the distribu- 
tion exponents tx corresponding to these ranges. 



III. CLUSTERS AND THEIR PROPERTIES 

We now proceed to apply the algorithms described in 
the previous section to the data presented in Table HI We 
start by providing some qualitative examples illustrating 
the complexity of the intermittent turbulent structures 
under study, as well as the performance of our clus- 
ter analysis code. Next, we present a detailed analysis 
of scaling behavior of these structures in the moderate- 
resolution runs I and III, followed by a comparative anal- 
ysis of the same properties in the high-resolution run II. 
Our primary objectives will be to identify relevant pa- 
rameters controlling the geometry of the observed struc- 
tures, to clarify the role of the initial conditions, and to 
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FIG. 3: {Color online) Geometric scaling of the current sheet 
structures detected in run I based on the j^ field, for sev- 
eral combinations of the snapshot time t and of the threshold 
parameter Oth'- i = 4, ath = 2 (medium gray or red) and 
ath ~ 3 (black); t — 5, ath = 2 (light gray or green) and 
ath ~ 3 (dark gray or blue). Dotted (dashed) vertical lines 
show the boundaries of the dissipative (inertial) scaling ranges 
used for computing the scaling exponents reported in Tables 
UlUlVI The statistics are roughly insensitive to the detection 
threshold ath and are stable in time. 



compare the geometry of the structures at inertial and 
dissipative scales. Finally, we discuss a possible link with 
self-organized criticality (SOC). 



A. The physical structures that emerge 

The top panel of Fig. [T] gives a perspective view of 
specific examples among about 700 dissipative regions 
detected for the lower resolution ABC flow (run I) at 
t = 4 and with ath = 2. The examples illustrate the 
complex multiscale nature of the j^ dissipation field, a 
typical feature of turbulent fluids. The current structures 
can be as large as the whole grid (not included in the 
figure to make smaller ones visible), or as small as several 
grid spacings. 

It should be emphasized that the upper panel in Fig. [1] 
is not produced by color-coding a continuous field as done 
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FIG. 4: {Color online) Scaling of probability distributions of 
current sheet structures detected in run I. Color coding and 
notations are the same as in the previous figure. Tilted dotted 
lines show power laws in the inertial and dissipative ranges 
for i = 4, ath = 2. As in the case of the geometric scaling, 
the shapes of the distributions are stable with respect to the 
detection threshold and time. 



often in turbulence visualizations. Each of the structures 
was first extracted by the algorithm described in the pre- 
ceding section. After that, roughly 1/10 of the structures 
were "placed back" in the domain according to their orig- 
inal positions and spatial orientations. We skipped the 
rest in order not to overcrowd the resulting picture. 

The apparent two-dimensional geometry of the struc- 
tures is typical for MHD turbulence, and it can be ob- 
served reliably over the entire inertial range of scales as 
we show in the next subsection. For smaller scales, the 
2D geometry becomes questionable, partly because the 
current sheets tend to fold or roll to form tubes which 
can no longer be resolved. 

The bottom panel of the same figure presents two typ- 
ical examples of large-scale dissipative regions extracted 
by our code, each occupying about 2000 grid nodes. As 
one can see, the regions may have rather complicated 
overall shapes associated with twisting and splitting of 
current sheets. Since the code does not rely on any a 
priory for the cluster shape or size, it can efficiently iden- 



TABLE II: Inertial range scaling exponents of current sheet 
structures in the runs with ABC and OT initial conditions 
(see Table |T]|; length scales for analysis are L £ [0.21, 1.30], 
t = 4, ath = 2, Af = 512. 



Dl 

TL 



Run I, j^ 



Run I, uj^ Run III, f 



1.00 ± 0.00 
2.20 ± 0.22 



1.00 ± 0.00 
2.57 ± 0.22 



1.00 ± 0.00 
1.86 ± 0.23 



1.00 ± 0.00 
1.94 ± 0.23 



Dr 0.94 ± 0.03 0.92 ± 0.01 0.93 ± 0.03 0.97 ± 0.02 
TR 2.14 ± 0.13 2.45 ± 0.24 1.90 ± 0.25 1.92 ± 0.21 



TABLE III: Dissipative range scaling exponents of current 
sheet structures in the runs with ABC and OT initial con- 
ditions (length scales L G [0.025,0.18], t = 4, ath = 2, 
AT = 512). 



Run I, j^ 



Run I, uj^ Run III, f Run III, lj^ 



Dl 1.00 ± 0.00 1.00 ± 0.00 1.00 ± 0.00 1.00 ± 0.00 

TL 0.75 ± 0.12 0.82 ± 0.04 1.31 ± 0.21 1.16 ± 0.08 

Dr 0.81 ± 0.04 0.84 ± 0.05 0.89 ± 0.02 0.82 ± 0.03 

TR 1.14 ± 0.09 1.22 ± 0.15 1.79 ± 0.09 1.65 ± 0.19 



Da 


1.71 ± 0.02 


1.72 ± 0.07 


1.72 ± 0.12 


1.92 ± 0.05 


Da 


1.40 ± 0.11 


1.44 ± 0.14 


1.53 ± 0.08 


1.29 ± 0.13 


TA 


1.66 ± 0.08 


1.81 ± 0.16 


1.66 ± 0.11 


1.39 ± 0.09 


TA 


0.95 ± 0.15 


1.04 ± 0.10 


1.38 ± 0.12 


1.29 ± 0.14 


Dv 


1.90 ± 0.06 


1.90 ± 0.03 


1.93 ± 0.12 


2.15 ± 0.05 


Dv 


1.51 ± 0.13 


1.54 ± 0.15 


1.60 ± 0.09 


1.40 ± 0.15 


TV 


1.61 ± 0.07 


1.61 ± 0.06 


1.49 ± 0.03 


1.36 ± 0.10 


TV 


1.10 ± 0.12 


1.16 ± 0.07 


1.39 ± 0.10 


1.40 ± 0.15 



Dh 0.18 ± 0.07 0.19 ± 0.04 0.21 ± 0.01 0.20 ± 0.05 
TH 3.51 ± 0.13 3.78 ± 0.09 4.13 ± 2.0 7.02 ± 2.0 



Dp 

Tp 



2.44 ± 0.10 
1.44 ± 0.06 



2.32 ± 0.04 
1.53 ± 0.05 



2.48 ± 0.17 
1.30 ± 0.08 



2.96 ± 0.08 
1.32 ± 0.06 



Dh 0.16 ± 0.01 0.15 ± 0.01 0.11 ± 0.02 0.11 ± 0.02 
TH 3.51 ± 0.13 3.78 ± 0.09 4.13 ± 2.0 7.02 ± 2.0 



Dp 

TP 



2.26 ± 0.14 
0.98 ± 0.08 



2.31 ± 0.13 
1.13 ± 0.08 



2.32 ± 0.08 
1.37 ± 0.08 



2.21 ± 0.24 
1.33 ± 0.08 



tify both simple and complex (e.g., folded or rolled-up) 
current sheets across the entire range of relevant spatial 
scales. 

Figure[2]further illustrates the ability of the code to de- 
tect complex dissipative structures. Here, the gray back- 
ground field represents the spatial distribution of j^ in x- 
y cross-sections of the data cube, with black (white) col- 
ors corresponding to the largest (smallest) current mag- 
nitude. The structures shown in the middle of the box 
(in red in the color version) correspond to a large-scale 
current sheet identified by our algorithm and embedded 
again into the turbulent flow to demonstrate its consis- 
tency with the surrounding MHD environment. In most 
of the slices, the highlighted current sheet consists of two 
disconnected pieces which only merge within a limited 
range of z values. Despite this topological complexity, 
the structure has been correctly identified as a single set 
of contiguous grid nodes. 



B. Statistics of structures 

Figure [3] shows the geometric scaling (dependence on 
length L, see Eq. ([S])) of the volume V, area A, Euclidian 
scale R, thickness H, dissipation rate Pj, and the com- 
plexity C on intermittent dissipative structures in the j^ 
field of run I. Figure U] shows the probability distribu- 
tions of the same parameters, except for C; the shade 
of gray (colors) used in the figures represent four differ- 
ent combinations of the time of the snapshots t (with 
t = 4 or t = 5, in units of the turn-over time of the 
problem, see [20I . |37|). and of the threshold ath (with 



ath = 2 or ath = 3). The results obtained for these 
parameters as well as for ath = 1 (not shown) are in- 
distinguishable within statistical uncertainty. Therefore, 
the scaling properties reported in this paper are not sen- 
sitive to the detection threshold, at least for ath G [1,3], 
and they do not vary on a time scale of the order of 
the turnover time. The dissipative and inertial ranges of 
scales are shown in both figures with vertical dotted and 
dashed lines respectively. 

The geometric scaling of the parameters V, A, R, and 
P exhibits clear-cut power-law behavior, with the log- 
log slopes (the D exponents of Eq. (fT5|)) undergoing 
slight changes at the transition between the two scal- 
ing regimes. The inertial values of Dy and Da suggest 
that the studied structures have a nearly two dimen- 
sional geometry. This is an expected result since the co- 
dimension of MHD turbulence is equal to unity, suggest- 
ing sheet-like structures embedded in three-dimensional 
space [32, l3j] ■ Similar Da estimates indicating sheet- like 
dissipative structures were obtained for the a; field in the 
same simulation run. A more careful inspection of the 
exponents obtained herein shows that the inertial range 
values of both Da and Dy are systematically below 2, 
suggesting that the structures have irregular edges. This 
interpretation is consistent with recently detected undu- 
lations of current sheet edges in the OT turbulence in 3D 

Unlike the geometric scaling of other size measures, the 
thickness H of the current sheets (Fig. [3]) does not vary 
significantly with L. It seems to saturate at the largest 
inertial scales revealing the existence of a characteristic 
thickness of the structures oi H = 0.04 — 0.05 (about 3- 



4 grid spacings). This thickness is hkely representative 
of the turbulent dissipation length (.diss whose estimated 
value is slightly larger than 3 x 27r/512 « 0.04. 

The geometric scaling of the topological complexity C 
also demonstrates a saturation in the inertial regime, and 
is clearly non-scale free at smaller scales. As already men- 
tioned in section II. C, the monotonic growth of C across 
a range of scales can be attributed to an increasingly 
complex shape of the structures, and is also affected by 
their folding. 

The distribution functions (Fig. |4]) demonstrate pro- 
nounced crossovers at the transition between the iner- 
tial and dissipative regimes. Overall, these crossovers arc 
more evident than the crossovers in the X{L) statistics 
shown in the previous figure. The thickness distribution 
is rather steep. It is likely to follow an exponential rather 
than a power-law decay, which is consistent with the exis- 
tence of a characteristic thickness H as discussed above. 

Tables |TT] and IIIII summarize inertial and dissipative 
range scaling exponents for runs I and III, using ath = 2, 
at i = 4, corresponding to the maximum of dissipation. 
The first column in each table refers to the log-log slopes 
of the red curves in Fig. |4] and Fig. [5l The exponent 
values reported in Table [ll] confirm that the geometry of 
both p and uP' structures observed in the inertial range is 
close to being two-dimensional. However, the dissipative 
range scaling (Table IIIip is significantly different, with 
the volume and area geometric exponents Dy and Da 
being close to 1.5, hinting at a fractal geometry of the 
structures with possible local anisotropy. 

The T exponents characterizing L (linear size) and R 
(Euclidian scale) arc almost identical in the inertial range 
(Table |ll| showing that cither parameter can be invoked 
as a measure of linear scale. On average, the identity 
L = R implies the absence of a preferred current or vor- 
ticity sheet orientation, indicating that at these scales, 
the MHD flows examined here are globally isotropic. The 
global isotropy obtains within the inertial range of scales 
and is not preserved at smaller scales, in accordance with 
previous results based on incompressible decaying MHD 
turbulence using ABC flows [SO]. Our analysis extends 
in an independent way these earlier flndings demonstrat- 
ing the inertial range global isotropy in both the OT and 
ABC runs. 

Note the values of th shown in Tables [Til and IIIII are 
equal. This is a result of the steep non-power law decay 
of the thickness distribution, which prevented us from 
measuring this exponent separately in the dissipative and 
inertial ranges of scales. 

When comparing the inertial distribution exponents 
characterizing the size of the structures, one can notice 
that these exponents arc systematically lower in run III 
(OT initial conditions) than in run I (ABC initial con- 
ditions). This difference is not large but is statistically 
significant for several exponent pairs. Alternatively, the 
dissipative range scaling shows the opposite effect (OT 
distribution exponents are higher than the ones in the 
ABC run). Since the dissipation range exponents are 



smaller on average than the inertial exponents, the iner- 
tial/dissipative distribution crossover is more pronounced 
in run I as compared to run III; on the other hand, this 
makes the statistics of OT structures closer to scale- 
invariant across the entire range of the scales, as evi- 
denced from the analysis of Dp, Dy, Tp, and Ty esti- 
mates shown in Tables |TT] and IIIII This may come from 
the fact that the OT flow has a well-defined structure 
with both partial zeros of the magnetic field (canceling 
of two components) and global zeroes (b = 0), leading 
to more ordered reconnection events and turbulence de- 
veloping at later times [S^l, whereas the ABC runs have 
some random noise added at small scales which leads to 
more wrinkled structures of lesser extent. 

To check the consistency of the pairs of D and t 
values we obtained, we also computed the exponents 
ax = Dx{tx — 1) which should be equal for all scale- 
invariant measures X of the examined sets of turbulent 
structures (due to the conservation of probability). We 
found the ax values to be approximately constant for 
most of the structure parameters, which confirms the va- 
lidity of our measurements. The only noticeable excep- 
tion is an whose value is inconsistent with the other a 
exponents. This can be seen as additional evidence of a 
non-power law scaling for the thickness of the dissipative 
structures, likely with a well-defined characteristic scale. 

Table HVl provides additional insight into the relation- 
ship between the various scaling exponents describing the 
parameters of j^ and w^ structures in the OT and ABC 
runs. The table shows estimated values of tlx for lin- 
ear size distribution exponents (based on the probability 
conservation) compared to the exponents tl evaluated 
directly from p{L) distributions. As one can notice, the 
values of tlx and r^ are in a reasonable agreement for 
the inertial range, but different by a roughly constant 
factor for the dissipative range. Also, the inertial range 
estimates tend to be lower for run III relative to run I; 
for the dissipative range they are higher. The difference 
between the exponents corresponding to the two scaling 
ranges seems to be less pronounced for the OT run. As an 
example, compare r^y = 2.2 (1.2) in the inertial (dissipa- 
tive) regimes of the ABC flow with the values ttv = 1.8 
(1.6) describing the same ranges in the OT flow. There- 
fore, the cross-over behavior in the linear size scaling is 
more evident for the ABC flow, in agreement with our 
conclusion based on the results in Tables [H] and IIIII 



This lack of complete universality in scaling of MHD 
flows can be related to similar findings in different con- 
texts. For example, it was shown in [40| that different 
power-law scaling for energy spectra can emerge with dif- 
ferent initial conditions for MHD flows having the same 
invariants {Et, He and Hm and Ey = Em at f = 0). 
Different energy spectra have also been observed in the 
presence of an imposed uniform and strong magnetic field 

ME!,!!!. 



TABLE IV: Comparison of linear size distribution exponents 
evaluated using the relation tlx = Dx{tx — 1) + i, X £ 
{R,A,V,P} (see Eqs. ([5|)- (|10[) ). with tl exponents obtained 
directly from p{L) distributions (f — 4, ath = 2, N — 512) for 
ABC (run I) and OT (run III). 



Run I, j Run I, 



Run III, f 



Inertial range 
TLR 2.07 

TLA 2.13 

TLV 2.16 

TLP 2.08 



{tlx) 2.11 ± 0.04 

TL 2.20 ± 0.22 

Dissipative range 

TLR 1.11 

TLA 0.92 

TLV 1.15 

TLP 0.96 

(tlx) 1.04 ± 0.11 

TL 0.75 ± 0.12 



2.34 
2.40 
2.15 
2.24 

2.28 ± 0.11 
2.57 ± 0.22 



1.19 
1.06 
1.25 
1.30 

1.20 ± 0.10 
0.82 ± 0.04 



1.84 
2.14 
1.94 
1.74 



1.92 
1.86 



± 0.17 
± 0.23 



1.70 
1.59 
1.63 
1.86 

1.70 ± 0.12 
1.31 ± 0.21 



1.89 
1.75 
1.78 
1.95 

1.84 ± 0.09 
1.94 ± 0.23 



1.54 
1.38 
1.56 
1.73 

1.55 ± 0.14 
1.16 ± 0.08 



C. Analysis at higher Reynolds number 

In order to determine what possible role the Reynolds 
number plays in the statistics of the structures studied 
in the preceding section, we now present an analysis of 
one snapshot for run II computed on a grid of 1536"^ 
points and taken at peak of dissipation, and contrast 
it with run I (the runs have Taylor Reynolds number 
of 1100 and 630, respectively). The high resolution run 
{N = 1536) is characterized by a larger Reynolds number 
and is therefore expected to generate more complex cur- 
rent and vorticity structures. A simple visual inspection 
of spatial patterns in j^ and w^ confirms this, and our 
quantitative analysis provides useful clues on the nature 
of the increased complexity in the high-resolution run. 

Figures [5] and [5] show comparative statistics of runs I 
and II for the j^ field. We expect, based on an approxi- 
mate convergence of j^ and w^ scaling exponents in the 
lower resolution runs (see Tables lllllIV|) . that the vor- 
ticity structures have a similar dependence on N. Note 
that the boundaries of inertial and dissipative ranges in- 
dicated by vertical lines in Figs. [5] and [6] are computed 
for N = 1536 and are thus different from the boundaries 
in the previous figures at lower Reynolds number. 

The geometric scaling (Fig. [5]) reveals a major dis- 
tinctive feature of run II: it has measurably thinner dis- 
sipative structures, as expected if we associate H with 
the dissipative scale idiss for this new run. Indeed, the 
volume of these structures is approximately half an or- 
der of magnitude smaller than the V estimates in run I 
made at the same linear scale L. At the same time, the 
scaling of the area, likely dependent on the integral scale 
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FIG. 5: {Color online) Geometric scaling of current structures 
in runs I and II detected by thresholding the j^ field at two 
different levels: A^ = 1536, ath = 2 (medium gray or red lines) 
and ath ~ 3 (black); A'^ = 512, ath = 2 (light gray or green) 
and ath ~ 3 (dark gray or blue). Dotted (dashed) vertical 
lines show the boundaries of the dissipative (inertial) scaling 
ranges. The higher Reynolds number ABC flow (run II) de- 
velops considerably thinner current sheet structures described 
by smaller volumes, and roughly the same surface areas com- 
pared to the lower resolution run. The dissipation rate in run 
II is slightly higher, and the geometry of current sheets gen- 
erated in this run is significantly more complex than the j^ 
structures observed in run I. 



of the flow Lint, is remarkably similar. The discrepancy 
between runs I and II has a straightforward explanation, 
namely, significantly thinner structures in the high reso- 
lution ABC run. On average, the values of H in this run 
are about three times smaller than the corresponding val- 
ues in run I (for the same L) . This difference is in agree- 
ment with the gain in the Reynolds number achieved due 
to the increase of the grid size from iV = 512 to 1536 (or 
the decrease of the viscosity, see Table |T| . As in the case 
of run I, the dimensionless thickness of the smallest in- 
ertial structures in run II is about 3-5 grid nodes. It is 
interesting that the scaling of the dissipation rate corre- 
sponding to structures for a higher Reynolds number flow 
is approximately the same as in run I, despite a signifi- 
cantly smaller volume and thickness of these structures. 
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FIG. 6: {Color online) Probability distributions of current 
structures in runs I and II. Tilted dotted lines show log-log 
regression slopes for run II, the remaining notations are the 
same as in the previous figure. With the exception of thick- 
ness distribution, all the statistics here exhibit power-law scal- 
ing with consistent sets of inertial range r exponents. 



Therefore, the current sheets generated on the grid with 
A^ = 1536"^ pomts are more intense (by a factor of « 3 
in terms of the current, usmg the H{L) scahng in the 
inertial regime). This is consistent with the finding that 
there is a finite dissipation rate in MHD in 3D [20| , a fact 
well-known in the 2D case |4j, |4J] and related to the pos- 
sibility of fast reconnection in MHD even in the absence 
of a Hall current (see |45| and references therein for recent 
developments) . Finally, the topological complexity of the 
structures in the high-resolution run is roughly twice the 
value characterizing the structures of run I, with C = 1 
matching the transition between the inertial and dissipa- 
tive regimes. A similar match was found for N = 512 
(see Fig. O. 

The probability distributions of j^ structures in runs I 
and II (Fig. ^ show that in spite of the essential differ- 
ence in the current sheet thickness, energy density and 
topological complexity evident from the geometric scal- 
ing, the probabilistic essence of the two runs is in fact 
quite similar. To a first approximation, the shape of the 
distributions is not influenced by Reynolds number, at 
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FIG. 7: (Color online) Comparison of scaling exponents of 
current structures in runs I and II (A'^ = 512 and A'^ = 1536, 
respectively) detected by thresholding the j^ field at the level 
of two standard deviations above the mean. The error bars 
are approximate confidence intervals (± three standard er- 
rors) for each exponent. The plots show consistent distribu- 
tion exponents but greater geometric exponents {Dl , Da and 
Dv) at higher Reynolds number. This difference becomes less 
prominent at higher detection thresholds (not shown). Simi- 
lar tendencies are observed for the vorticity structures. 



least for the small range examined here. As in the case 
of the lower Reynolds number simulation at lower res- 
olution, the distributions in run II exhibit well-defined 
"breaks" coinciding with the transition between the two 
scaling regimes (and again shifted toward smaller scales 
for the high Ry run) . It also shows a rapid decay of the 
probability distribution p{H) , consistent with the equiv- 
alent distribution in run I. 

Figure [7] displays the numerical values of several key 
scaling exponents, reflecting the differences and the sim- 
ilarities between high- and low-Reynolds number ABC 
runs. The overlapping error bars indicate that the mean 
exponents are indistinguishable at the 99% confldence 
level. The plots confirm that the distribution exponents 
are roughly independent of Ry whereas the geometric 
exponents tend to be larger in run II and are also closer 
to the value of two, as expected for idealized purely two- 
dimensional structures. The difference is consistent with 
the previously discussed observation that the current 
sheets in run II are considerably thinner and might there- 
fore be better described by a 2D scaling model within the 
inertial range of scales. 

Finally, note that the proximity of D^ in run II to a 
value of unity indirectly indicates that the current sheets 
generated in this fluid are somewhat more isotropic in 
terms of their global (large-scale) spatial orientation com- 
pared to run I. For a fully isotropic macroscopic oricnta- 
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FIG. 8: [Color online) Scaling characteristics of current and 
vorticity structures in ABC and OT flows (runs I and III, 
A*' = 512, i = 4, ath = 2) addressing tfie role of the velocity 
- magnetic field alignment. In run I, structures in j'^ and 
oj^ are shown with light gray (green) and dark gray (blue) 
lines, whereas in run III, they are given with middle gray 
(red) and black lines. Note a similarity in the volume and 
area statistics, and a significant difference in the alignment 
patterns of the two flows (x is the cosine of the local alignment 
angle between v and b). In addition to a stronger normalized 
He, the OT flow exhibits systematically thinner structures 
(smaller thickness for a given linear size, for both p and uj^) 
compared to the ABC flow. 



tion, we expect the two linear measures, R and L, to be 
directly proportional. 

Overall, we conclude that the Reynolds number of the 
flow, at least in the case of unit magnetic Prandtl number 
(v = r]), influences the geometry of the resulting dissipa- 
tive structures, making them more 2D-like (higher aspect 
ratio) and better mixed in space, but it does not seem to 
alter their statistical properties as reflected by the family 
of the distribution functions we have studied. 



D. The role of velocity-magnetic field correlations 

It has been known for quite some time that the amount 
of correlation between the velocity and the magnetic held 
plays a significant role in the dynamics of MHD turbu- 
lence (see \4M for a review, and more recent works in 



[23 . |47|). This role can be global, altering the scaling of 
energy spectra when He, normalized by Et, is strong, 
i.e., close to ±1; it can also be important, even when 
globally weak, since structures with strong alignment be- 
tween V and b develop rapidly in a turbulent flow |2l| 
(note that He is not globally positive definite). It thus 
appears as a natural application of our detection algo- 
rithm to examine the properties of the selected clusters 
in this light. 



Our examination of the data reveals a significant dif- 
ference in the velocity - magnetic field alignment for the 
OT and ABC runs, in agreement with earlier analyses 
[2l| . Defining x as the cosine of the local alignment an- 
gle between v and b, the OT run is characterized by an 
asymmetric p(x) distribution (with a skewness of « 1.0), 
having a sharp maximum at x = +1- The alignment 
distribution suggests a prevailing parallel orientation of 
V and b fields inside the inertial range structures, with 
an average x ~ 0.4. The ABC run shows a nearly sym- 
metric p{x) distribution (skewness « 0) with the average 
alignment close to zero as well (see Fig. |8l bottom right 
panel). The stronger velocity - magnetic field alignment 
in the OT fluid may be the primary reason for its dis- 
tinct scaling behavior as represented in Tables |ll] - IIVI 
As we have already stressed, the OT r exponents tend to 
be less sensitive to the crossover between the inertial and 
dissipative ranges of turbulence compared to the ABC 
run. This tendency is especially clear in the statistics 
of vorticity structures (black lines in Fig. [51 see also the 
last column in Table HI)) . At the same time, the geometric 
scaling of most of the parameters oi j^ and w^ structures, 
in particular V{L) and A{L) dependencies, is practically 
the same for the two flows except for the geometric scal- 
ing of H, implying that the current and vorticity sheets 
are somewhat thinner in the OT run. 



Indeed, the alignment effect appears to have limited 
or no impact on scaling behavior of the intermittent 
structures. Surprisingly, the differences in the geomet- 
ric and probabilistic scaling of OT and ABC turbulence 
are strikingly small compared to the dramatic difference 
in the x distributions characterizing the two runs. It is 
also interesting that the p{x) distributions have the same 
functional form for current and vorticity structures (for 
a given initial condition), suggesting the existence of a 
strong i - Lo coupling correlated with the v - 6 align- 
ment; this result may be linked to the fact that, when 
writing the MHD equations in terms of the Elsasser fields 
z = V ± b for which the nonlinear terms reduce to an 
advection of one field by the other, one sees that the dy- 
namics strongly couples the velocity and magnetic fields 
(and their derivatives) [22|. Other alignments could be 
considered [43| , from the point of view of structure anal- 
ysis, basically those having a direct impact on the dy- 
namics, such as the Lamb vector v x w, the Lorentz force 
j X b and Ohm's law v x b. This is left for future work. 
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IV. A POSSIBLE LINK WITH 
SELF-ORGANIZED CRITICALITY 

The analysis presented in this section is motivated by 
the intensively debated connection between intermittent 
structures in turbulence and SOC, the latter having been 
discussed intensely in the literature over the last decade. 
In particular, SOC has been proposed as the underlying 
physical mechanism responsible for the intermittency of 
the dissipation field in high-Reynolds number turbulent 
fluids [50r - t53i | (for a review of SOC in the context of Solar 
Wind and the magnetosphere, see, e.g., [5J,[53]). It has 
been suggested that dissipativc regions can communicate 
over large distances, by analogy with critical avalanches 
in sandpile models of SOC, producing conditions for a 
statistically steady state of nonequilibrium critical be- 
havior resp onsible for multifractal inhomogeneous dissi- 
pation [56|-|58|. 

To test the SOC avalanche hypothesis, one needs to 
obtain a collection of probability distribution functions 
describing the dissipative regions, and to study their 
scaling. The hallmark of SOC is the power-law shape 
of avalanche distributions over a number of parameters, 
some of which are studied here. According to the defi- 
nition of SOC, avalanches are essentially spatiotemporal 
objects composed of all grid nodes involved in the for- 
mation of a given dissipative structure over its entire life 
cycle. Consequently, in order to rigorously verify SOC 
behavior in a turbulent fluid, one needs to detect and 
analyze its dissipative structures in both space and time. 
The results of our present analysis refer to static three- 
dimensional vorticity and current clusters observed at a 
fixed time; thus, these results cannot provide a definite 
answer to the question of whether or not incompressible 
MHD turbulence is related to SOC. Nevertheless, they 
can be used to make some preliminary estimates (see |6l| 
for an analysis of flaring activity in one space dimension, 
and [62[ for the 2D case). 

In the following derivation, we are assuming that the 
spatial intermittent structures explored in the previous 
sections are static snapshots of dynamic intermittent 
events evolving in space and time. Using the SOC ap- 
proach, each of these events can be described by the 
spatiotemporal size 5" representing the total number of 
grid nodes involved in the event over its life time T. 
The distributions of S and T are expected to scale as 
p{S) ^ S'-^s and p{T) - T~^^ . The avalanche size and 
lifetime scaling exponents {ts and tt) are usually con- 
sidered to be the primary measures of criticality defined 
by the universality class of a particular set of symmetries 
describing local interactions between the nodes [63| . 

Due to the absence of temporal dimension in our 
present analysis, neither of the two SOC exponents is 
directly accessible. However, by applying once again the 
conservation of probability, we can evaluate them indi- 
rectly through 



TABLE V: Avalanche size and avalanche lifetime scaling ex- 
ponents estimated in the dissipative range using the relations 
rsx = l + Dx{tx - l)/Ds and ttx = 1 + Dx{tx - 1)/Dt, 
with X G {A,V,P} (see Eq. ^E^) with the mean-field geo- 
metric exponents Ds ~ 2 and Dt ~ 1 and the tx exponents 
taken from Table III (t = 4, atn = 2, iV = 512) 





Run I, f 


Run I, uj'^ 


Run III, f 


Run III, uj'^ 


TSA 


0.97 


1.03 


1.29 


1.19 


rsv 


1.08 


1.12 


1.31 


1.28 


TSP 


0.98 


1.15 


1.43 


1.36 


(rs) 


1.01 ± 0.06 


1.10 ± 0.06 


1.34 ± 0.08 


1.28 ± 0.09 


tta 


0.93 


1.06 


1.58 


1.37 


ttv 


1.15 


1.25 


1.62 


1.56 


ttp 


0.96 


1.30 


1.86 


1.73 


{tt) 


1.01 ± 0.12 


1.20 ± 0.13 


1.69 ± 0.15 


1.55 ± 0.18 



in which S ^ L^^, T ^ L^^ , and X stands for one of 
the static measures of the structures exhibiting power- 
law scaling as already discussed. 

We start by analyzing in this light the dissipative sub- 
range of scales. In the case of the ABC flow, the prox- 
imity of Tx to the value 1 in that range (see Table IIII[) 
makes uncertainty of Ds and Dt unimportant. For a 
wide range of Ds and Dt estimates presented in the 
SOC literature, and for various choices of X, Eq. p^ 
predicts that the value of both the avalanche size and 
the lifetime distribution exponents for this flow are also 
close to unity. Thus, for instance, by plugging in Ty and 
Dy for current structures in the ABC flow [64|, and us- 
ing mean-field values for Ds and Dt (respectively 2 and 
1) [65[, one gets ts ~ tt ~ 1- Interestingly enough, 



TS,T = 1 + Dx{tx -l)/Ds^T, 



(16) 



the same calculation for the OT run yields a significantly 
different result, namely rg « 1.3 and tt ~ 1.6. 

Table IVl summarizes the estimated values of ts and tt 
exponents obtained using the relation (|16p in which we 
plug in the dissipative range scaling exponents character- 
izing volume, surface area, and energy dissipation rate in 
j^ and uj^ structures. We do not use the linear size expo- 
nents in this calculation since they are less reliable due 
to their dependence on the orientation of the structures 
(sec the discussion of the small-scale anisotropy in section 
III.B). 

While the ABC avalanche exponents obtained from 
Eq. p6p in the dissipative range arc somewhat low com- 
pared to SOC exponents usually reported in the liter- 
ature, the OT exponents clearly fall within the range 
of values expected for many SOC sandpiles. Thus, for 
example, they are a very close match to the 2D real- 
ization of the directed Abelian sandpile model (DASM) 
[6a | , an exactly solvable version of the paradigmatic Bak- 
Tang-Wiesenfeld model [ll]. The DASM avalanche dis- 
tributions arc described by the exponents Ts — 4/3 and 
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ory of critical phenomena, the tendency seen in Fig. [9] is 
iisuaUy approximated by truncated power-law distribu- 
tions of the form 



pix) = X--- fix/x,) , 



(17) 



10 10 10 10 10 10 10 
Volume 



FIG. 9: (Color online) Probability distributions of the di- 
mensionless volume and area in runs I (middle gray or red), 
II (light gray or green), and III (dark gray or blue), demon- 
strating a wider range of small-scale power-law behavior in 
the high Rv (high-resolution) simulation. 



tt = 3/2 which are nearly the same as the OT values re- 
ported in Table V. Several other avalanching models are 
approximately consistent with the predicted OT expo- 
nents, such as, e.g., the Manna two-state model [63, the 
Bak-Snappen model of punctuated evolution |68ll. and 
the absorbing state phase transition model |70|, all in 
two spatial dimensions. The distinctive feature of DASM 
which might be responsible for the best match with criti- 
cal behavior in the dissipative range of the OT flow is the 
existence of the preferred spatial direction in which the 
avalanche fronts propagate. As we have mentioned ear- 
lier in the text, the local spatial anisotropy seems to be 
a significant factor in the dissipative range scaling of the 
studied fluids, and therefore the observed agreement be- 
tween OT and DASM avalanche exponents may be more 
than just a coincidence. 

Besides the existence of well-known SOC universality 
classes consistent with the exponents obtained from the 
analysis, a significant piece of evidence for the involve- 
ment of SOC dynamics in the formation of multifractal 
intermittent dissipation fields comes from the statistics 
of dimensionless measures of dissipative structures rep- 
resenting their size and geometry in terms of discrete grid 
nodes, by analogy with sandpile simulations. A predicted 
effect for a SOC system is an expansion of the range of 
the power-law scaling of such measures for increasingly 
large lattice sizes, known as the critical finite-size scaling 
(FSS) behavior. 

Our analysis shows that the probability distributions 
of dimensionless quantities of turbulent structures (e.g., 
dimensionless volume V/6^ and area A/ 6"^; see Fig. [9]) ex- 
hibit signatures of FSS, assuming that SOC avalanches 
do develop in the dissipative (sub-inertial) range of scales. 
A comparison of dimensionless distributions in high- and 
low-resolution runs shows that the range of power-law 
scaling describing the smallest structures, and thus in- 
volving a limited number of grid nodes, expands towards 
larger scales as N increases from 512 to 1536. This type 
of behavior is among the most distinctive signatures of 
SOC systems. It indicates that the intrinsic mechanisms 
of avalanching dynamics are scale- free with no limit other 
than the limited size of the Reynolds number. In the the- 



where / is an appropriate scaling function and Xc the 
(apparent) characteristic scale of structures increasing 
with a discrete system size N as Xc ^ iV"^. Due to 
the scale-free nature of SOC states, we also expect that 
Kx is close to the corresponding geometric exponent Dx 
(because the largest linear scale associated with A'' has 
the same effect on the distributions of structure sizes as 
the intermediate scales associated with L). However, it 
should be recognized that unlike FSS in ordinary SOC 
simulations with simple boundary conditions, the upper 
scale of the presumed SOC behavior in our turbulent 
runs is controlled by a complex process - the inertial 
range turbulent cascade. Consequently, the functional 
form of the scaling function / in a turbulent fluid should 
perhaps include a scale-free component accounting for 
the fluid turbulence at larger scales (in contrast, e.g., to 
the exponentially decaying / commonly used in sandpile 
simulations). 

The interpretation of SOC exponents found in the in- 
ertial range is more ambiguous. Technically, they are 
not far from the values ts ~ 1.5 and tt = 2 describ- 
ing SOC dynamics in 3-dimensional stochastic and de- 
terministic directed sandpiles (see [69| and refs. therein). 
However, this similarity can be misleading as the above- 
mentioned models possess a distinct geometry consist- 
ing of two spatial dimensions in which dissipative events 
can grow isotropically and one dimension allowing for 
unidirectional (directed) growth only. Whether or not 
the growth of dissipative structures in the inertial-range 
MHD turbulence contains such a preferred direction im- 
plying strong mesoscopic anisotropy, remains to be veri- 
fied. Until then, we only associate the dissipative range 
with SOC behavior. It may be the case that the accuracy 
of the present analysis in the inertial range is insufficient, 
or that SOC behavior is only limited to the dissipative 
range. A direct spatio-temporal analysis or growing and 
decaying dissipative structures will be instrumental for 
validating our SOC observations in the inertial range and 
for reducing the uncertainties in the exponents. Such 
analysis will also allow for a study of whether the ergod- 
icity assumptions often used in the study of turbulent 
flows are valid |83j . 

Overall, the results obtained in this paper suggest 
that small-scale dissipative structures, observed below 
the smallest inertial range scale, are associated with SOC. 
If this hypothesis is correct, the small-scale intermittency 
in 3D MHD turbulence can be interpreted as a propaga- 
tion of local instabilities from small to large scales indica- 
tive of SOC avalanches. This propagation should reflect a 
tendency of the smallest dissipative structures, such as, 
e.g., current filaments, to merge into larger clusters in 
an avalanche fashion, and it does not necessarily imply a 
transport of energy in Fourier space in the opposite direc- 
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tion as to the main (direct) turbulent cascade of energy. 
One could think for example of larger dissipation events, 
like major flares in the solar corona, emerging from a 
cooperative behavior of smaller events (e.g., nanoflares), 
an effect reminiscent of non-locality of nonlinear interac- 
tions between widely separated scales [SOj. This possi- 
bility has been discussed intensively in the literature, in 
particular in the framework of forest-fire models of mul- 
tifractal inhomogeneous dissipation in turbulent media 
[7l| . Our results seem to be the first (but so far indi- 
rect) evidence of such behavior in incompressible MHD 
observed through DNS in 3D, even though, of course, 
SOC behavior in MHD has been advocated for a long 
time followin g th e pioneering paper of Lu and Hamilton 
[ll (seealso |5l[5l|63,[73l). 

In a more general context, to fully describe SOC behav- 
ior in turbulence one would have to understand several 
fundamental aspects shaping the dynamics of the inter- 
mittent structures. For example, what plays the role of 
a threshold in MHD turbulence, so that avalanches (dis- 
sipative events) of various sizes can happen? It may be 
a collective effect triggered by sweeping of large scales, 
pushing together magnetic field lines of opposite polar- 
ities to come in close contact as has been proposed by 
Klimas et al. [T^l, or it may be the instability of cur- 
rent (or vorticity) sheets below a given thickness (at a 
fixed viscosity /resistivity). If SOC is indeed identified, 
other relevant questions involve, e.g., what is the un- 
derlying phase transition, what is the order parameter, 
and what field can be associated with the susceptibility 
that diverges at the critical point? These questions may 
not have clear answers, as the identification of critical- 
ity in fiuid and MHD turbulence is not straightforward, 
although an example in the context of atmospheric pre- 
cipitation in tropical convection has been put forward 
recently [49| . 

We can speculate on one of the possible sources 
of burstiness in the dissipation and its relation with 
avalanches. Starting from the pioneering work of On- 
sager, Lee and Kraichnan |74| - [76| . one can think of the 
dynamical evolution of a turbulent flow being due to non- 
linear interactions with weak forcing and weak dissipa- 
tion balancing each other. Solutions of the ideal trun- 
cated equations obtain at late times in the simplest in- 
stance equipartition between all the modes, with zero en- 
ergy flux. At intermediate times and intermediate scales, 
one observes turbulent dynamics with non-zero flux [77[ , 
the "dissipation" of large-scale energy being associated to 
a turbulent eddy viscosity due to the thermalized modes 
at small scale (see jTSl . uM fo'^ similar results for helical 
flows, and for 2D MHD). This dynamics have also been 
observed in viscous cases, e.g., in Navier-Stokes at high 
resolution, where the resulting flow can be decomposed 
into a set of coherent structures with a spectrum close 
to Kolmogorov, and a large number of modes at small 
scales in thermal equilibrium [80| . 

Related with these results, it has been known for some 
time, and in different instantiations of turbulent flows. 



that the energy flux of a given sign on average, has in fact 
huge fluctuations of both signs and of amplitudes much 
larger than the mean (of order unity for characteristic 
velocities and lengths [/q = 1 and Lq = 1, respectively; 
see for example ^] and references therein, and [S^l for 
studies of regions with zero flux in models of turbulence) . 
These large fluctuations in the flux can be attributed to 
the balance between forcing and dissipation mentioned 
above, and to the two components (one thermalized and 
random, one turbulent and coherent) identified in turbu- 
lent flows at small scales. The interplay between the two 
components can result in a bursty flux transfer of energy 
to the small scales, as observed in particular when look- 
ing at dissipation and reconnection events [ij]. These 
bursts are the needed excursions that lead the system 
away from equilibrium and may give rise to a state of 
criticality, in order to dissipate the energy accumulated 
over various lapses of time through the injection mecha- 
nism. Some of these events will trigger other events, by 
pushing around structures that through their associated 
pressure fields can make contact with other structures 
that may in turn destabilize. 

Finally, it is worth mentioning that intermittent bursts 
of energy transfer have been observed in MHD in the 
so-called shell models for turbulent flows [8J]. These 
models can be viewed as a poor-man template, set on 
a lattice, for the temporal evolution of the Navier-Stokes 
or MHD equations (see [8g and references therein for 
a recent review); in the simplest case, one retains only 
nearest-neighbor interactions which are built in such a 
way that the quadratic invariants are preserved. Shell 
models have been known to exhibit avalanche behavior 
as well (see [8a,|86[ for a discussion), although SOC inter- 
pretation of energy avalanches in such simplified models 
of turbulent cascade remains questionable, as there are 
indications that only spatially-distributed systems with 
clear time-scale separation between the driving and dis- 
sipation mechanisms are able to exhibit robust SOC sig- 
natures ISTl. 



V. CONCLUSION 
A. Summary of the results 

In this paper, we have analyzed three sets of data 
stemming from high-Reynolds number numerical sim- 
ulations of MHD turbulence in three dimensions at a 
magnetic Prandtl number of unity; periodic boundary 
conditions are assumed and there is no imposed uni- 
form magnetic field nor is there any forcing. Initial con- 
ditions are either the so-called 3D Orszag-Tang vortex 
(OT) or the Arn'old-Befframi-Childress (ABC) fiow; the 
two flows have different velocity-magnetic field correla- 
tion He- Numerical resolutions range from 512'^ grid 
points to 1536'^, with Taylor Reynolds numbers R\ vary- 
ing from 630 to 1100. 

We find that current and vorticity sheets behave in 
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similar fashion and that, overaU, the probabihstic prop- 
erties of these structures do not depend on either Rx or 
the normahzed value of He, i.e., the degree to which v 
and b are globally aligned. Other factors like the detec- 
tion threshold used for the analysis (within 1 to 3 stan- 
dard deviations above the mean) or the time at which 
the turbulent fluid is analyzed (after reaching the peak 
of the dissipation) appear to be irrelevant as well. 

As expected, dissipative structures are thinner and 
more complex (and closer to being two-dimensional) the 
higher the Reynolds number. Compared to the ABC 
flow, the OT run has more efficient kinetic energy dis- 
sipation (higher dissipation rate per unit volume) and 
its distributions of structure parameters are closer to a 
"mo no fractal" shape across all scales, i.e., approximately 
described by a single power-law. The high Rx run has 
a higher dissipation rate per unit volume (but roughly 
the same rate per unit surface) . and essentially the same 
form of probability distributions as in the low Rx run. 

The incrtial scaling exponents characterizing the OT 
and ABC flows are similar, with differences between the 
exponents more marked at smaller scales. The exponents 
obtained at these scales (in the dissipative subrange) can 
be associated with SOC universality classes. Our findings 
suggest that whereas the inertial-range scaling is more 
likely to be dictated by MHD turbulent cascades (energy 
spectra, and structure functions in general), the dissipa- 
tive range of scales may be governed by self-organized 
criticality; this is perhaps the main finding that stems 
from our cluster analysis. 

Why is the SOC scaling found in the dissipation range? 
Simply because this is where the approximately ideal dy- 
namics breaks down; at the small-scale end of the inertial 
range detailed conservation of quadratic invariants is vi- 
olated by non-zero dissipation. In fact, if there are singu- 
larities in such flows, the dissipation can be order unity 
(in terms of the characteristic velocity and length). Such 
exchanges must be bursty insofar as they are concen- 
trated on a small scale; they are rare as they only occur 
in special cases, and thus they must be strong as they 
provide the dissipation on average to balance the energy 
injected by the larger scales. In other words, the stochas- 
tic (and irreversible) element comes from the fact that 
dissipative events occur where vortices and currents of 
opposite polarities meet, at random time because of the 
randomness of large-scale structures for times longer than 
the eddy-turnover time (see, e.g., [83 for studies of 1// 
noise in turbulent flows). This results in strong bursts 
when the system undergoes dissipation events (such as 
reconnection) with reconfiguration of the fields. In this 
simplified picture, the critical parameter for transition 
may be the Kolmogorov dissipation length Idiss, in scale- 
space, or rather the dimensionlcss local Reynolds number 
uii/v of order unity, by definition, at € ~ idiss- 



B. Final remarks 



In conclusion, we have shown that the cluster algo- 
rithm presented in section II. B and used throughout this 
study can easily and systematically detect a multitude 
of current and vorticity structures in a turbulent flow, 
with a large dynamical range in their intensity, and that 
it readily leads to an analysis of their relevant physical 
properties. As important perhaps, one can then study 
some other statistics of such structures, such as the align- 
ment between velocity and magnetic field. Although the 
analysis presented in this paper deals with properties of 
decaying fiows in the vicinity of the peak of dissipation, 
one can note that such flows are thought to behave sim- 
ilarly to the statistically steady forced case, since their 
dynamics is quasi-steady for some interval of time around 
that peak (sec, e.g., [SO])- It would however be of interest 
to study the statistically stationary case of MHD turbu- 
lence; this is left for future work. 

Relating the scaling exponents found in a given analy- 
sis to other, more traditional measures of complexity in 
a turbulent flow dealing, e.g., with correlation functions, 
is not necessarily a straightforward task. Some fascinat- 
ing results concerning the behavior of two-dimensional 
Navier-Stokes turbulence have been unraveled recently 
jGOj with no clear direct connection to the scaling of, say, 
the energy spectrum; in this 2D case with an inverse cas- 
cade of energy to large scales, the study of the zero-line 
vorticity contours led to the discovery of a link with a spe- 
ciflc class of percolation and anomalous diffusion through 
the scaling laws for, e.g., length versus diameter. The fact 
that the dynamics of turbulent flows contains elements of 
critical phenomena and conformal invariance associated 
with invariance properties and symmetry groups of the 
underlying equations points to the need to further our 
studies of such flows using scaling tools. 

There are other algorithms that can examine coher- 
ent structures in turbulent flows, following the pioneering 
work for 2D Navier-Stokes fluids [24| . Prominent among 
them nowadays is the wavelet decomposition [89| which 
leads to a vision of turbulent flows in three dimensions as 
a set of coherent structures (vortex filaments in the fluid 
case) with a Kolmogorov spectrum, together with inco- 
herent quasi- Gaussian eddies which contain most of the 
degrees of freedom and which are in some sense slaved to 
the coherent structures; a similar analysis has been done 
recently in MHD [80| . The use of visualization techniques 
and lossless compression of data is also a possible tool to 
analyze turbulent structures |34l . |90[. Combining such 
tools with the analysis of hypercubes of data taking into 
account the temporal dimension of structures may prove 
fruitful, but in three dimensions this represents a chal- 
lenge that we want to tackle in the near future, both 
for fluids and MHD; it will allow for a better connec- 
tion between turbulence, intermittency, structures and 
self-organized criticality. 
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